#Declarations

var N[y.end.sim, n.class], A[n.class,n.class,y.end.sim], rN[y.end.sim,n.class], f.neg[y.end.sim], f.pos[y.end.sim], f.con[y.end.sim], N.total[y.end.sim], phi[y.end.sim], r.sero.calf[y.end.sim], r.sero.yrlg[y.end.sim], r.sero.adult[y.end.sim]
#Data assingnment for initial conditions
data {
	mu.1 <-y.N[1]
	
	}
model{

#Priors 
#<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-
#Coeficients in recruitment regresion
#from analysis of YNP data



#Priors from Fuller paper
b[1] ~ dbeta(77,18.08)
b[2] ~ dbeta(36.5,20.54)
b[3] ~ dbeta(3.2, 11.36)


#Priors on pi and rho for serology test, from meta analysis
pi ~ dbeta(476,42)
rho ~ dbeta(403,15467)


#shape parameters for gamma distribution of census standard deviations

alpha.var ~ dgamma(.001,.001)
nu.var ~ dgamma(.001,.001)

#prior on sighting probability
p.sight ~ dbeta(370,9.50)


#prior on vaccination efficacy
eta~ dbeta(4,4)  #efficacy of vaccine

#Survival
#Juvenile
p[1] ~ dunif(0,1)
#Adult. Mean and sd of survival taken from YNP white binder
#mean<-.96, s[d <- .03, converted to shape parameters for beta
p[2] ~ dbeta(40, 1.667) 
p[3] ~ dunif(0,1)

#weakly informative prior on offspring sex ratio (mean <- .5, sd <- .05)
m ~ dbeta(49.5,49.5)

#disease parameters

beta ~ dunif(0,50)

v ~ dbeta(y.shapes.v[1],y.shapes.v[2])

psi ~ dbeta(y.shapes.psi[1], y.shapes.psi[2])


#probability of seropositive in removals
for(t in 1:y.end.sim){
	r.sero.calf[t] ~ dbeta(1,1)
	r.sero.yrlg[t] ~ dbeta(1,1)
	r.sero.adult[t] ~ dbeta(1,1)
	}

#stochasticity in process model

	sigma.p ~ dunif(0,5)
	tau.p <- 1/(sigma.p^2)
	# tau.p ~ dgamma(.01,.01)
	# sigma.p <- 1/sqrt(tau.p)
	

#End of priors
#<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-

###inital conditions
#Initial conditions
for(j in 1:n.class){
	shape.d[j] <- 1
}
comp.1 ~ ddirch(shape.d[1:n.class])
tau.1 <- 1/var.N.obs
N.1 ~ dnorm(mu.1,tau.1)
for(j in 1:n.class){
	N.start[j]<- N.1*comp.1[j] 
}	

N[1,1:n.class] <- N.start[1:n.class]


#inital conditions for dervied quantities and latent states
# total population size

f.neg[1] <- b[1]
f.pos[1] <- b[2]
f.con[1]  <- b[3]

#z.calf[1] ~ dunif(0,100)  #These are needed to allow initial values for t > 1
#z.yrlg[1] ~ dunif(0,100)
#z.adult[1] ~ dunif(0,100)

diff.pos <- f.neg[1]-f.pos[1]
diff.con <- f.neg[1]-f.con[1]
N.total[1] <-sum(N[1,1:n.class])

#prevalance by age class
	pv.cow[1] <-  sum(N[1,6:7])/(sum(N[1,6:7])+N[1,3])
	pv.yrl[1] <- N[1,5] /(N[1,2] + N[1,5])
	pv.calf[1] <- N[1,4]/(N[1,4]+N[1,1]) 
#Population compostion 

	calf.ratio[1] <- (N[1,1] + N[1,4]) / N.total[1] 
	cow.ratio[1] <- (N[1,3] + sum(N[1,6:7]))/ N.total[1] 
	yrl.ratio[1] <- (N[1,2] + N[1,5])/ N.total[1]
	bull.ratio[1] <- N[1,8]/N.total[1]
	
#Probability of transmission

#Mixedt transmission, z is random variable
z~dbeta(1,1)
#initial value of probability of transmission to allow monitoring, not used for inference
phi[1] <- 1-exp(-beta * (p[2]*N[1,6] + p[2]*psi*N[1,7]) / ((1-m)*(p[1]*N[1,1]+ p[1]*N[1,4]) + p[2]*(sum(N[1,2:3]) + p[2]*sum(N[1,5:7])))^z)
#proportion of adult cows that are infectious
in.per[1] <- (N[1,6]/(sum(N[1,6:7]) + N[1,3]))
#propotion of sero-postive adult cows that are infectious
in.per.pos[1] <- (N[1,6])/(sum(N[1,6:7]))


Re[1] <-0 #placeholder for first value in vector

#initialize new infection count and proportions

#initialize new infection count and proportions

NewI[1,1] <-0
NewI[1,2] <- 0
NewI[1,3] <-0
NewI[1,4] <-0
NewI[1,5] <-0
NewI[1,6] <- 0
NewI[1,7] <- 0
NewI[1,8] <- 0
NewIP[1,1] <-0
NewIP[1,2] <- 0
NewIP[1,3] <-0
NewIP[1,4] <-0
NewIP[1,5] <-0
NewIP[1,6] <- 0
NewIP[1,7] <- 0
NewIP[1,8] <- 0
NewIP[1,9] <- 0
#proportion of new infections by vertical transmission

v.b[1] <- (1-m)*p[2]* v *(f.con[1]*N[1,6] + f.pos[1]*N[1,7])
vert.per[1] <- v.b[1] / sum(N[1,4:6])


#estimate proportion of each age class in removals that is seropositive
for(i in 1:length(y.remove.index)){
	
	
	z.calf[i] ~ dbin(r.sero.calf[y.remove.index[i]],y.rsero.calf[y.remove.index[i],2])
	y.rsero.calf[y.remove.index[i],1] ~ dbin((pi*z.calf[i] + rho*(y.rsero.calf[y.remove.index[i],2] - z.calf[i]))/y.rsero.calf[y.remove.index[i],2],y.rsero.calf[y.remove.index[i],2])
#		
	
	#z.yrlg[y.remove.index[i]] ~ dbin(r.sero.yrlg[y.remove.index[i]],y.rsero.yrlg[y.remove.index[i],2])
	#y.rsero.yrlg[y.remove.index[i],1] ~ dpois(pi*z.yrlg[t] + rho*(y.rsero.yrlg[t,2] - z.yrlg[t]))
	z.yrlg[i] ~ dbin(r.sero.yrlg[y.remove.index[i]],y.rsero.yrlg[y.remove.index[i],2])
 	y.rsero.yrlg[y.remove.index[i],1] ~ dbin((pi*z.yrlg[i] + rho*(y.rsero.yrlg[y.remove.index[i],2] - z.yrlg[i]))/y.rsero.yrlg[y.remove.index[i],2],y.rsero.yrlg[y.remove.index[i],2])
#	
	# z.adult[t] ~ dbin(r.sero.adult[t],y.rsero.adult[t,2])
	# y.rsero.adult[t,1] ~ dpois(pi*z.adult[t] + rho*(y.rsero.adult[t,2] - z.adult[t]))
	
	z.adult[i] ~ dbin(r.sero.adult[y.remove.index[i]],y.rsero.adult[y.remove.index[i],2])
	y.rsero.adult[y.remove.index[i],1] ~ dbin((pi*z.adult[i] + rho*(y.rsero.adult[y.remove.index[i],2] - z.adult[i]))/y.rsero.adult[y.remove.index[i],2],y.rsero.adult[y.remove.index[i],2])
}

#initial conditions for the suvival vector--no removals occcured at time 1 

s[1,1] <- p[1]	
s[1,2] <- p[2]	
s[1,3] <- p[2]	
s[1,5] <- p[2]
s[1,6] <- p[2]
s[1,7] <- p[2]
s[1,8] <- p[3]




#constant for portion of the year elapsed before peak of removals
	
q<-3/4


#Process Model<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-

for(t in 2:y.end.sim){

##cacluate recrutiment rates for sero positive and sero negative females
#Recruitment parameters
f.neg[t] <- b[1]
f.pos[t] <- b[2]
f.con[t]  <- b[3]




#removals
#####Matrix version of removals
#estimate age / sex composition of removals
	r.comp[t,1:4]~ddirch(y.rage[t,1:4])
	#r.sero.calf[t] ~ dbeta(y.rsero.calf[t,1], y.rsero.calf[t,2])
	#r.sero.yrlg[t] ~ dbeta(y.rsero.yrlg[t,1], y.rsero.yrlg[t,2])
	#r.sero.adult[t] ~ dbeta(y.rsero.adult[t,1], y.rsero.adult[t,2])
	
	#probabilities that a removed animal is seropositive	
		
	#seronegtives removed
	#calves
	rN[t,1] <-y.removals[t]*r.comp[t,1]*(1-r.sero.calf[t])
	#yearlings
	rN[t,2] <-y.removals[t]*r.comp[t,2]*(1-r.sero.yrlg[t])
	#adults
	rN[t,3] <-y.removals[t]*r.comp[t,3]*(1-r.sero.adult[t])
	#######sero positives removed
	#calves
	rN[t,4] <-y.removals[t]*r.comp[t,1]*r.sero.calf[t]
	#yearlings
	rN[t,5] <-y.removals[t]*r.comp[t,2]*r.sero.yrlg[t]
	#removals of seropositive adult cows must be partitioned among 6, and 7
	rN[t,6] <-y.removals[t]*r.comp[t,3]*r.sero.adult[t] * N[t-1,6]/sum(N[t-1,6:7]) 	
	rN[t,7] <-y.removals[t]*r.comp[t,3]*r.sero.adult[t] * N[t-1,7]/sum(N[t-1,6:7])  	
	
	
	
	#bulls removed
	rN[t,8] <- y.removals[t]*r.comp[t,4]
	
#suvival probability based on removals and natural mortality. 
	s[t,1] <- max(0, p[1]-p[1]^(1-q)*rN[t,1]/N[t-1,1])
	s[t,2] <- max(0, p[2]-p[2]^(1-q)*rN[t,2]/N[t-1,2])
	s[t,3] <- max(0, p[2]-p[2]^(1-q)*rN[t,3]/N[t-1,3])
	s[t,4] <- max(0, p[1]-p[1]^(1-q)*rN[t,4]/N[t-1,4])
	s[t,5] <- max(0, p[2]-p[2]^(1-q)*rN[t,5]/N[t-1,5])
	s[t,6] <- max(0, p[2]-p[2]^(1-q)*rN[t,6]/N[t-1,6])
	s[t,7] <- max(0, p[2]-p[2]^(1-q)*rN[t,7]/N[t-1,7])
	s[t,8] <- max(0, p[3]-p[3]^(1-q)*rN[t,8]/N[t-1,8])



##frequency depdendent transmission
numer_phi[t]<-s[t,6]*N[t-1,6] + s[t,7]*psi*N[t-1,7]
denom_phi[t] <- ((1-m)*(s[t,1]*N[t-1,1]+s[t,4]*N[t-1,4]) + sum(s[t,2:3]*N[t-1,2:3]) + sum(s[t,5:7]*N[t-1,5:7]))^z

phi[t] <- 1-exp(-beta * numer_phi[t] /denom_phi[t])


#####Transition Matrix
#row 1
A[1,3,t]<-s[t,3]*f.neg[t]*(1-phi[t])
A[1,6,t]<-s[t,6]*f.con[t]*(1-v)*(1-phi[t])
A[1,7,t]<-s[t,7]* (f.pos[t] *(1-psi)*(1-phi[t])  +  f.pos[t] * psi*(1-v)*(1-phi[t]))


#row 2
A[2,1,t]<-s[t,1]*(1-m)*(1-phi[t])
#row 3
A[3,2,t]<-s[t,2]*(1-phi[t])
A[3,3,t]<-s[t,3]*(1-phi[t])

#row4
A[4,3,t]<-s[t,3]*f.neg[t]*phi[t]
A[4,6,t]<-s[t,6]*f.con[t]*(v+phi[t]-v*phi[t])
A[4,7,t]<-s[t,7]*(f.pos[t]*psi*(v+phi[t]-v*phi[t]) + f.pos[t]*(1-psi)*phi[t])

#row5
A[5,1,t]<-s[t,1]*phi[t]*(1-m)
A[5,4,t]<-s[t,4]*(1-m)

#row6
A[6,2,t]<-s[t,2]*phi[t]
A[6,3,t]<-s[t,3]*phi[t]
A[6,5,t]<-s[t,5]
A[6,7,t] <-0

#row 7
A[7,6,t]<-s[t,6]
A[7,7,t]<-s[t,7]

#row 8
A[8,1,t]<-s[t,1]*m
A[8,4,t]<-s[t,4]*m
A[8,8,t]<-s[t,8]



#estimates of median of posterior distribution
mu[t,1:n.class] <- A[,,t]%*%N[t-1,1:n.class]

#get logs of deterministic predictions
for(j in 1:8){
	log_mu[t,j] <- log(max(1,mu[t,j]))
	}

log_N[t,1:8] ~ dmnorm(log_mu[t,1:8],tau.p*y.I) 

#exponentiate log of population size
for(j in 1:8){
	N[t,j] <- min(exp(log_N[t,j]),10000)
	}


N.total[t] <- sum(N[t,1:n.class])


#Population compostion 

	calf.ratio[t] <- (N[t,1] + N[t,4]) / N.total[t] #calves
	cow.ratio[t] <- (N[t,3] + sum(N[t,6:7])) / N.total[t]
 #adult females
 	yrl.ratio[t] <- (N[t,2] + N[t,5]) / N.total[t]

	bull.ratio[t]<-N[t,8]/N.total[t]
	
	

#sero-prevalace by age
	pv.cow[t] <-  sum(N[t,6:7])/(sum(N[t,6:7])+N[t,3])
	pv.yrl[t] <- N[t,5] /(N[t,2] + N[t,5])
	pv.calf[t] <- N[t,4]/(N[t,4]+N[t,1]) 
	
#proportion of adult cows that are infectious
in.per[t] <- (N[t,6]/(sum(N[t,6:7]) + N[t,3]))
#propotion of sero-postive adult cows that are infectious
in.per.pos[t] <- (N[t,6])/(sum(N[t,6:7]))
#proportion of new infections by vertical transmission

v.b[t] <- s[t,2]* v * f.con[t]*N[t-1,6] + s[t,7]*f.pos[1]*N[t-1,7]
vert.per[t] <- v.b[t] / sum(N[t,4:6])




##Effective reproductive ratio
#Re[t] <- ((1-m) *N[t,4] + phi[t] * p[2]*sum(N[t,1:3])) / (N[t-1,6]+psi*N[t-1,7])



#Effective reproductive ratio
Re[t] <- ((1-m) * N[t,4] + phi[t]*sum(s[t-1,1:3]*N[t-1,1:3])) / (N[t-1,6]+psi*N[t-1,7])
#Effective rerproductive ratio in the absence of managment removals
Re2[t] <- ((1-m) * N[t,4] + phi[t]*(p[1]*N[t-1,1]+p[2]*(N[t-1,2]+N[t-1,3]))) / (N[t-1,6]+psi*N[t-1,7])



#number of new infections

NewI[t,1] <- (1-m)*N[t-1,6]*s[t,6]*f.con[t]*v # vertical from from infectious
NewI[t,2] <- (1-m)*N[t-1,7]*s[t,7]*psi*f.pos[t]*v # vertical from recrudesced



NewI[t,3] <- (1-m)*(N[t-1,3]*s[t,3]*f.neg[t] + N[t-1,7]*s[t,7]*f.pos[t]) * phi[t] #horizontal juvenile
NewI[t,4] <- N[t-1,2]*s[t,2]*phi[t] # yearling horizontal 
NewI[t,5] <- N[t-1,3]*s[t,3]*phi[t] # adult horizontal
NewI[t,6] <- N[t-1,7]*s[t,7]*psi #recrudesence
NewI[t,7] <- sum(NewI[t,1:2]) #total vertical
NewI[t,8] <- sum(NewI[t,3:5]) #total horizontal

#proportions of new infections
NewIP[t,1:6] <- NewI[t,1:6] / sum(NewI[t,1:6]) #1:6 because this is the total number of new infections
NewIP[t,7] <-sum(NewI[t,1:2])/sum(NewI[t,1:6]) #proportion vertical
NewIP[t,8] <-sum(NewI[t,3:5])/sum(NewI[t,1:6]) #proportion horizontal, remainder is proportion recrudesence.
NewIP[t,9] <- NewI[t,6]/sum(NewI[t,1:6]) #proportion recrudensed

} #end process model


#averages across years
xcalf.ratio <- mean(calf.ratio[6:41])
xcow.ratio <- mean(cow.ratio[6:41])
xyrl.ratio <- mean(yrl.ratio[6:41])
xbull.ratio <- mean(bull.ratio[6:41])

xpv.cow <- mean(pv.cow[6:41])
xpv.yrl <- mean(pv.yrl[6:41])
xpv.calf <- mean(pv.calf[6:41])

xin.per <- mean(in.per[6:41])
xin.per.pos <- mean(in.per.pos[6:41])
xvert.per <- mean(vert.per[6:41])

xphi <- mean(phi[6:41])

for(j in 1:8){
	xNewIP[j]<-mean(NewIP[6:41,j])
}


#Data models<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-

#Likelihood for annual counts################################
##Estimate of observation error in census data

for(i in 1:length(y.N.varindex[])){
	#estimate variance based on observed variance
	y.N.var[y.N.varindex[i]] ~ dgamma(alpha.var,nu.var)
}

#Estimate of the mean observation variance

var.N.obs <- (alpha.var/nu.var)

for(t in 1:length(y.N.index)){
	#method of moments to get gamma shape parameters
	#from model estimates
	alpha[t]<- N.total[y.N.index[t]]^2/var.N.obs
	nu[t] <- N.total[t]/var.N.obs
	lambda[t] ~ dgamma(alpha[t],nu[t]) 
	lambda.cor[t]<-p.sight*lambda[t]
	y.N[t] ~ dpois(lambda.cor[t])
	
	#simulated data for posterior predictive check
	y.N.new[t] ~ dpois(p.sight*lambda[t])
	y.N.sq[t] <- (p.sight*lambda[t] - y.N[t])^2
	y.N.new.sq[t] <- (p.sight*lambda[t]-y.N.new[t])^2
}
	#Bayesian P value
	P.N.mean <- step(mean(y.N.new[]) - mean(y.N[]))
	P.N.sd <- step(sd(y.N.new[]) - sd(y.N[]))
	P.N.sq <- step(sum(y.N.new.sq[]) - sum( y.N.sq[]))
	
##  Likelihoods for population composition data ##############################################


#proportion of calves in population from aerial counts
for(t in 1:length(y.iratio.calf)){    #
	mu.calf[t] <- (N[y.iratio.calf[t],1] + N[y.iratio.calf[t],4]) / N.total[y.iratio.calf[t]]
#calcuate shape parameters for beta from data and estimates from external, hierachical model of ratio.calf data, y.sigma.calf[t] and ratio.calf[t] which are read in as data
	shape1.calf[t] <- max(.0001,(mu.calf[t]^2-mu.calf[t]^3-mu.calf[t]*y.sigma.calf[t]^2)/y.sigma.calf[t]^2)
	shape2.calf[t] <- max(.0001, (mu.calf[t]-2*mu.calf[t]^2+mu.calf[t]^3-y.sigma.calf[t]^2+mu.calf[t]*y.sigma.calf[t]^2)/y.sigma.calf[t]^2)
	y.ratio.calf[t] ~ dbeta(shape1.calf[t], shape2.calf[t])
	#simulate new data for posterior predictive checks
	y.ratio.calf.new[t] ~ dbeta(shape1.calf[t], shape2.calf[t])
	y.ratio.calf.sq[t] <- (mu.calf[t] - y.ratio.calf[t])^2
	y.ratio.calf.new.sq[t] <- (mu.calf[t] - y.ratio.calf.new[t])^2

} 
#Bayesian P values

P.ratio.calf.mean <- step(mean(y.ratio.calf.new[])-mean(y.ratio.calf[]))
	P.ratio.calf.sd <- step(sd(y.ratio.calf.new[])-sd(y.ratio.calf[]))
	P.ratio.calf.sq <- step(sum(y.ratio.calf.new.sq[]) - sum(y.ratio.calf.sq[]))


#proportions of yearling cows, adult cows, and non-yearling bulls in the population from ground counts
#for the first two years, the yearlings were lumped with the adult females. 
for(t in 1:2){
	alpha.p[t,1] <- bull.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,2] <- (cow.ratio[y.ground.index[t]] + yrl.ratio[y.ground.index[t]])*y.alpha_0[t]
	alpha.p[t,3] <- 0.0001*y.alpha_0[t]
	alpha.p[t,4] <- calf.ratio[y.ground.index[t]]*y.alpha_0[t]
	y.p.mu[t,1:4] ~ ddirch(alpha.p[t,1:4])
	}
for(t in 3:length(y.ground.index[])){
	alpha.p[t,1] <- bull.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,2] <- cow.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,3] <- yrl.ratio[y.ground.index[t]]*y.alpha_0[t]
	alpha.p[t,4] <- calf.ratio[y.ground.index[t]]*y.alpha_0[t]
	y.p.mu[t,1:4] ~ ddirch(alpha.p[t,1:4])
	}






#Likelihoods for serology data####################################

#serology of calves:
for(j in 1:length(y.isero.calf[])){
		
	y.pos.calf[j] ~ dbin((N[y.isero.calf[j],4]*pi + N[y.isero.calf[j],1]*rho)/ (N[y.isero.calf[j],4]+N[y.isero.calf[j],1]), y.n.calf[j])
	
	#simulate data for posterior predictive check
	y.pos.calf.new[j] ~ dbin((N[y.isero.calf[j],4]*pi + N[y.isero.calf[j],1]*rho)/ (N[y.isero.calf[j],4]+N[y.isero.calf[j],1]), y.n.calf[j])

	y.pos.calf.sq[j] <- (y.pos.calf[j]-pv.calf[j])^2
	y.pos.calf.sq.new[j] <- (y.pos.calf.new[j] - pv.calf[j])^2
			
	}
	#Bayesian P value
	P.sero.mean[1] <- step(mean(y.pos.calf.new[])-mean(y.pos.calf[]))
	P.sero.sd[1] <- step(sd(y.pos.calf.new[])-sd(y.pos.calf[]))
	P.sero.sq[1] <- step(sum(y.pos.calf.sq.new[])-sum(y.pos.calf.sq[]))

	
#serology of yearlings
for(j in 1:length(y.isero.yrl[])){
	#sero-prevalance in calves
	
	y.pos.yrl[j] ~ dbin((N[y.isero.yrl[j],5]*pi + N[y.isero.yrl[j],2]*rho)/ (N[y.isero.yrl[j],5]+N[y.isero.yrl[j],2]), y.n.yrl[j])
	
	#simulate data for posterior predictive check
	y.pos.yrl.new[j]~ dbin((N[y.isero.yrl[j],5]*pi + N[y.isero.yrl[j],2]*rho)/ (N[y.isero.yrl[j],5]+N[y.isero.yrl[j],2]), y.n.yrl[j])

	y.pos.yrl.sq[j] <- (y.pos.yrl[j]-pv.yrl[j])^2
	y.pos.yrl.sq.new[j] <- (y.pos.yrl.new[j] - pv.yrl[j])^2
	}
	
	#Bayesian P value
	P.sero.mean[2] <- step(mean(y.pos.yrl.new[])-mean(y.pos.yrl[]))
	P.sero.sd[2] <- step(sd(y.pos.yrl.new[])-sd(y.pos.yrl[]))
	P.sero.sq[2] <- step(sum(y.pos.yrl.sq.new[])-sum(y.pos.yrl.sq[]))

	
#serology of adult cows
for(j in 1:length(y.isero.cow[])){
	y.pos.cow[j] ~ dbin((sum(N[y.isero.cow[j],6:7])*pi + N[y.isero.cow[j],3]*rho)/ (sum(N[y.isero.cow[j],6:7])+N[y.isero.cow[j],3]), y.n.cow[j])

	
	#Simulate data for posterior predictive check
	y.pos.cow.new[j] ~ dbin((sum(N[y.isero.cow[j],6:7])*pi + N[y.isero.cow[j],3]*rho)/ (sum(N[y.isero.cow[j],6:7])+N[y.isero.cow[j],3]), y.n.cow[j])
	y.pos.cow.sq[j] <- (y.pos.cow[j]-pv.cow[j])^2
	y.pos.cow.sq.new[j] <- (y.pos.cow.new[j] - pv.cow[j])^2

	}
	#Bayesian P value
	P.sero.mean[3] <- step(mean(y.pos.cow.new[])-mean(y.pos.cow[]))
	P.sero.sd[3] <- step(sd(y.pos.cow.new[])-sd(y.pos.cow[]))
	P.sero.sq[3] <- step(sum(y.pos.cow.sq.new[])-sum(y.pos.cow.sq[]))


#prediction of phi
#based on momment matching for beta distribution:

for(t in 1:length(y.iphi)){
	mu.phi[t] <- phi[y.iphi[t]]
	sq.err1[t] <- (y.phi[t] - mu.phi[t])^2
	sq.err2[t] <- (( y.phi[t] - mu.phi[t])^2)/y.phi.sd[t]^2
#	#expressions below were througly tested to assure they produce propper shape parameters from mean and sd's
	a.phi[t] <-max(.001,(mu.phi[t]^2-mu.phi[t]^3-mu.phi[t]*y.phi.sd[t]^2)/y.phi.sd[t]^2)
	b.phi[t] <- max(.00001,(mu.phi[t]-2*mu.phi[t]^2+mu.phi[t]^3-y.phi.sd[t]^2+mu.phi[t]*y.phi.sd[t]^2)/y.phi.sd[t]^2)
	Py.phi[t] <- dbeta(y.phi[t], a.phi[t],b.phi[t])
}

mspe1 <- sum(sq.err1)/length(y.iphi)
mspe2 <- sum(sq.err2)/length(y.iphi)
P_yphi_phi <- prod(Py.phi[])

}
/*End of Model*/